Depth From Disparity

CSC 471 Computer Vision

Pooneet Thaper and Jin Zhang

Professor Zhigang Zhu

Attempt 1: Living Room

Overview

Images were taken of my living room each aligned with the edge of a table and offset by a baseline of 25cm.

Materials

• A small tripod and a phone camera mount adapter were used to provide stability
• A tape measure was used to measure 25 centimeter intervals from the left edge of the table
• A phone (Oneplus 6) was used to capture the photos with default settings

Procedure

• Five intervals were measured (0cm, 25cm, 50cm, 75cm, 100cm from the left edge of the table)
• The left edge of the tripod was aligned with each of these intervals
• The front edge was aligned with the front edge of the table for each photo
• Images were mistankenly taken with autofocus which could have caused different focal lengths for each

Discussion

• The focus (found afterwards by looking at the metadata of the photos) was found to be 4.25mm on all images
• Disparity equation can be used to find the depth of each area of the image: $$ Z = D = f \frac{B}{d} $$ where: $$ D = \textrm{depth} $$ $$ f = \textrm {focal length} $$ $$ B = \textrm{baseline} $$ $$ d = \textrm{disparity} = x_r - x_l $$

import os
from PIL import Image
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

if on_colab:
  image_dir = 'Depth-from-Disparity/images/living_room/'
else:
  image_dir = '../images/living_room/'

images_available = sorted(os.listdir(image_dir), key= lambda image_path: int(image_path.split("_")[-2]))

images = [Image.open(os.path.join(image_dir, image_path)).rotate(180) for image_path in images_available]
images = [image.resize((800,600)) for image in images]

# All 5 images
num_images = len(images)
fig, ax = plt.subplots(nrows=1, ncols=num_images, sharey=True, figsize=(20, 3))

for axis, image, image_path in zip(ax, images, images_available):
    axis.imshow(image)
    axis.set_title("Baseline {} cm".format(image_path.split("_")[-2]))

plt.show()

# The first and last image with a baseline of 100cm (1m)
fig, ax = plt.subplots(nrows=1, ncols=2, sharey=True, figsize=(20, 7))

for axis, image, image_path in zip(ax, [images[0]] + [images[-1]], [images_available[0]] + [images_available[-1]]):
    axis.imshow(image)
    axis.set_title("Baseline {} cm".format(image_path.split("_")[-2]))

plt.show()

# The second and fourth image with a baseline of 50cm
chosen_images = [1, 3]

fig, ax = plt.subplots(nrows=1, ncols=2, sharey=True, figsize=(20, 7))

for axis, image, image_path in zip(ax, [images[chosen_images[0]]] + [images[chosen_images[1]]],
                                   [images_available[chosen_images[0]]] + [images_available[chosen_images[1]]]):
    axis.imshow(image)
    axis.set_title("Baseline {} cm".format(image_path.split("_")[-2]))

plt.show()

chosen_grayscale_images = [get_intensity_image(images[i])/255 for i in chosen_images]

fig, ax = plt.subplots(nrows=1, ncols=2, sharey=True, figsize=(20, 7))

for axis, image, image_path in zip(ax, chosen_grayscale_images, 
                                   [images_available[chosen_images[0]]] + [images_available[chosen_images[1]]]):
    axis.imshow(image, cmap='gray')
    axis.set_title("Normalized Grayscale Baseline {} cm".format(image_path.split("_")[-2]))

plt.show()

# Checking to see if the correspondence is working
i = 200
j = 150

kernel_Y_min = i - margin
kernel_Y_max = i + margin + 1
kernel_X_min = j - margin
kernel_X_max = j + margin + 1
print(kernel_Y_min,kernel_Y_max,kernel_X_min,kernel_X_max)

min_Y_window = max(i - search_range, 0)
max_Y_window = min(i + search_range + 1, chosen_grayscale_images[0].shape[0])
min_X_window = j - margin
max_X_window = min(j + (2*search_range) + 1, chosen_grayscale_images[0].shape[1])
print(min_Y_window, max_Y_window, min_X_window, max_X_window)

192 209 142 159
175 226 142 201

# Area from right image to be matched in left image
plt.imshow(chosen_grayscale_images[1][kernel_Y_min:kernel_Y_max, kernel_X_min:kernel_X_max], cmap='gray')
plt.title('Area to be matched from right image')
plt.show()

# Area in left image to be searched for right image portion
fig,ax = plt.subplots(ncols=2, nrows=1, figsize=(16, 7))

ax[0].imshow(chosen_grayscale_images[0][min_Y_window:max_Y_window, min_X_window:max_X_window], cmap='gray')
ax[0].set_title('Area to be searched in left image')

ax[1].imshow(chosen_grayscale_images[1][min_Y_window:max_Y_window, min_X_window:max_X_window], cmap='gray')
ax[1].set_title('Same area in the right image (to show disparity)')
plt.show()

# Output of correspondence (second coordinate is discrepency)
output = correspondence_2D(chosen_grayscale_images[1][kernel_Y_min:kernel_Y_max, kernel_X_min:kernel_X_max], 
                               chosen_grayscale_images[0][min_Y_window:max_Y_window, min_X_window:max_X_window])
max_coords = np.unravel_index(np.argmin(output, axis=None), output.shape)
print(max_coords)

(17, 8)

# Check the returned area against the original
fig,ax = plt.subplots(ncols=2, nrows=1, figsize=(16, 7))

ax[0].imshow(chosen_grayscale_images[1][kernel_Y_min:kernel_Y_max, kernel_X_min:kernel_X_max], cmap='gray')
ax[0].set_title('Original region in right image')
ax[1].imshow(chosen_grayscale_images[0][min_Y_window+max_coords[0]:min_Y_window+max_coords[0]+(kernel_Y_max-kernel_Y_min), 
                                      min_X_window+max_coords[1]:min_X_window+max_coords[1]+(kernel_X_max-kernel_X_min)], cmap='gray')
ax[1].set_title('Corresponding region in left image')
plt.show()

Attempt 2: City College

Overview

Images were taken of City College from across the street each aligned with the edge of a window and offset by a baseline of 11 inches (the length of a peice of paper).

Materials

• A window was used to provide stability and alightment between the photos
• A piece of computer paper was used to measure 11 inch intervals from the left edge of the window
• A phone (Oneplus 6) was used to capture the photos with default settings (having a focal length of 4.25mm)

Procedure

• Five intervals were measured (0in, 11in, 22in, 33in, and 44in from the left edge of the window)
• The bottom edge of the phone was aligned with each of these intervals
• The phone was kept flat against the glass of the window to keep the optical axes parallel

Discussion

• The focus (found afterwards by looking at the metadata of the photos) was 4.25mm on all images
• Disparity equation can be used to find the depth of each area of the image: $$ Z = D = f \frac{B}{d} $$ where: $$ D = \textrm{depth} $$ $$ f = \textrm {focal length} $$ $$ B = \textrm{baseline} $$ $$ d = \textrm{disparity} = x_r - x_l $$

import os
from PIL import Image
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

if on_colab:
  image_dir = 'Depth-from-Disparity/images/city_college/'
else:
  image_dir = '../images/city_college/'

images_available = sorted(os.listdir(image_dir), key= lambda image_path: int(image_path.split("_")[-2]))

images = [Image.open(os.path.join(image_dir, image_path)).rotate(180) for image_path in images_available]
images = [image.resize((800,600), resample=Image.LANCZOS) for image in images]

# All 5 images
num_images = len(images)
fig, ax = plt.subplots(nrows=1, ncols=num_images, sharey=True, figsize=(20, 3))

for axis, image, image_path in zip(ax, images, images_available):
    axis.imshow(image)
    axis.set_title("Baseline {} cm".format(image_path.split("_")[-2]))

plt.show()

# The first and last image with a baseline of 44in
fig, ax = plt.subplots(nrows=1, ncols=2, sharey=True, figsize=(20, 7))

for axis, image, image_path in zip(ax, [images[0]] + [images[-1]], [images_available[0]] + [images_available[-1]]):
    axis.imshow(image)
    axis.set_title("Baseline {} cm".format(image_path.split("_")[-2]))

plt.show()

import numpy as np

def correspondence_2D(array_a, array_b, center_y = 0, center_x = 0):
  # convolves array_a across array_b but subtracts rather than multiplies
  # assumes dimensions of array_b greater than or equal to corresponding of array_a

  results = np.zeros(array_b.shape[:2])
  convs = 0

  for y in range(len(array_b)):
    for x in range(len(array_b[0])):
      if convs % 100 == 0:
        pass
      convs += 1
      for i in range(len(array_a)):
        for j in range(len(array_a[0])):
          try:
            try:
              a = array_a[i,j]
            except:
              a = 0
            try:
              idx_y = y+i-center_y
              idx_x = x+j-center_x
              if idx_y < 0 or idx_x < 0:
                b = 0
              else:
                b = array_b[idx_y,idx_x]
            except:
              b = 0
            results[y,x] += np.absolute(a-b)
            #print(y,x,i,j,y+i-center_y,x+j-center_x,a,b)
          except:
            pass
  return results

def correspondence_3D(array_a, array_b):
  # convolves array_a across array_b but subtracts rather than multiplies
  # assumes dimensions of array_b greater than or equal to corresponding of array_a
  
  results = np.zeros(array_b.shape[:2])
  convs = 0
  
  for y in range(len(array_b)):
    for x in range(len(array_b[0])):
      if convs % 100 == 0:
        pass
      convs += 1
      for i in range(len(array_a)):
        for j in range(len(array_a[0])):
          for k in range(3):
            try:
              try:
                a = array_a[i,j,k]
              except:
                a = 0
              try:
                b = array_b[y+i,x+j,k]
              except:
                b = 0
              results[y,x] += np.absolute(a-b)
              
            except:
              pass
  return results

def NTSC_Luminance_Equation(pixel):
    luminance_coefficients = [0.299, 0.587, 0.114]
    return np.dot(pixel, luminance_coefficients)

def get_intensity_image(image):
    img_array = np.array(image)
    return np.array([[NTSC_Luminance_Equation(img_array[i, j, :]) 
                    for j in range(img_array.shape[1])] 
                   for i in range(img_array.shape[0])])

# The second and third image with a baseline of 11in seem to be good candidates
chosen_images = [1, 2]

fig, ax = plt.subplots(nrows=1, ncols=2, sharey=True, figsize=(20, 7))

for axis, image, image_path in zip(ax, [images[chosen_images[0]]] + [images[chosen_images[1]]],
                                   [images_available[chosen_images[0]]] + [images_available[chosen_images[1]]]):
    axis.imshow(image)
    axis.set_title("Baseline {} cm".format(image_path.split("_")[-2]))

plt.show()

chosen_grayscale_images = [get_intensity_image(images[i])/255 for i in chosen_images]

fig, ax = plt.subplots(nrows=1, ncols=2, sharey=True, figsize=(20, 7))

for axis, image, image_path in zip(ax, chosen_grayscale_images, 
                                   [images_available[chosen_images[0]]] + [images_available[chosen_images[1]]]):
    axis.imshow(image, cmap='gray')
    axis.set_title("Normalized Grayscale Baseline {} cm".format(image_path.split("_")[-2]))

plt.show()

kernel_size = 17 #assumed odd since that is usual
margin = kernel_size//2
search_range = 25
log_file = open('log.txt','a+')

# Checking to see if the correspondence is working
i = 400
j = 450

kernel_Y_min = i - margin
kernel_Y_max = i + margin + 1
kernel_X_min = j - margin
kernel_X_max = j + margin + 1
print(kernel_Y_min,kernel_Y_max,kernel_X_min,kernel_X_max)

min_Y_window = max(i - search_range, 0)
max_Y_window = min(i + search_range + 1, chosen_grayscale_images[0].shape[0])
min_X_window = j - margin
max_X_window = min(j + (2*search_range) + 1, chosen_grayscale_images[0].shape[1])
print(min_Y_window, max_Y_window, min_X_window, max_X_window)

392 409 442 459
375 426 442 501

# Area from right image to be matched in left image
plt.imshow(chosen_grayscale_images[1][kernel_Y_min:kernel_Y_max, kernel_X_min:kernel_X_max], cmap='gray')
plt.title('Area to be matched from right image')
plt.show()

# Area in left image to be searched for right image portion
fig,ax = plt.subplots(ncols=2, nrows=1, figsize=(16, 7))

ax[0].imshow(chosen_grayscale_images[0][min_Y_window:max_Y_window, min_X_window:max_X_window], cmap='gray')
ax[0].set_title('Area to be searched in left image')

ax[1].imshow(chosen_grayscale_images[1][min_Y_window:max_Y_window, min_X_window:max_X_window], cmap='gray')
ax[1].set_title('Same area in the right image (to show disparity)')
plt.show()

# Output of correspondence (second coordinate is discrepency)
output = correspondence_2D(chosen_grayscale_images[1][kernel_Y_min:kernel_Y_max, kernel_X_min:kernel_X_max], 
                               chosen_grayscale_images[0][min_Y_window:max_Y_window, min_X_window:max_X_window])
max_coords = np.unravel_index(np.argmin(output, axis=None), output.shape)
print(max_coords)

(19, 29)

# Check the returned area against the original
fig,ax = plt.subplots(ncols=2, nrows=1, figsize=(16, 7))

ax[0].imshow(chosen_grayscale_images[1][kernel_Y_min:kernel_Y_max, kernel_X_min:kernel_X_max], cmap='gray')
ax[0].set_title('Original region in right image')
ax[1].imshow(chosen_grayscale_images[0][min_Y_window+max_coords[0]:min_Y_window+max_coords[0]+(kernel_Y_max-kernel_Y_min), 
                                      min_X_window+max_coords[1]:min_X_window+max_coords[1]+(kernel_X_max-kernel_X_min)], cmap='gray')
ax[1].set_title('Corresponding region in left image')
plt.show()

import csv
i = margin
j = margin
step = 5

log_written = []
log = []
current_row = 0
current_col = 0

while i < chosen_grayscale_images[1].shape[0]-margin:
    while j < chosen_grayscale_images[1].shape[1]-margin:
        if i != current_row or (current_col-j) >= 200:
            print(i,j)
            current_row = i
            current_col = j

        kernel_Y_min = i - margin
        kernel_Y_max = i + margin + 1
        kernel_X_min = j - margin
        kernel_X_max = j + margin + 1
        #print(kernel_Y_min,kernel_Y_max,kernel_X_min,kernel_X_max)

        min_Y_window = max(i - search_range, 0)
        max_Y_window = min(i + search_range + 1, chosen_grayscale_images[0].shape[0])
        min_X_window = j - margin
        max_X_window = min(j + (2*search_range) + 1, chosen_grayscale_images[0].shape[1])
        #print(min_Y_window, max_Y_window, min_X_window, max_X_window)

        output = correspondence_2D(chosen_grayscale_images[1][kernel_Y_min:kernel_Y_max, kernel_X_min:kernel_X_max], 
                                   chosen_grayscale_images[0][min_Y_window:max_Y_window, min_X_window:max_X_window])
        max_coords = np.unravel_index(np.argmin(output, axis=None), output.shape)
        #print(max_coords)
        log.append((i,j,
                    kernel_Y_min,kernel_Y_max,kernel_X_min,kernel_X_max,
                    min_Y_window, max_Y_window, min_X_window, max_X_window,
                    max_coords[0], max_coords[1]))
        #print()
        j += step
        
    if len(log):
        with open('log.csv','a',newline='') as file:
        writer = csv.writer(file)
        for entry in log:
            writer.writerow(entry)
        log_written += log
        log = []
        file.close()
        
    i += step
    j = margin

148 413
153 8
158 8
163 8
168 8
173 8
178 8
183 8
188 8
193 8
198 8
203 8
208 8
213 8
218 8
223 8
228 8
233 8
238 8
243 8
248 8
253 8

parallel = Image.open('../images/parallel_programming.png')
plt.figure(figsize=(20,10))
plt.imshow(parallel)
plt.show()

import csv
log_file = open('test/log.csv','r')
reader = csv.reader(log_file)
discrepancies = np.array([[np.NaN]*chosen_grayscale_images[0].shape[1]]*chosen_grayscale_images[0].shape[0])
depths = np.array([[np.NaN]*chosen_grayscale_images[0].shape[1]]*chosen_grayscale_images[0].shape[0])

points = []
values = []

for row in reader:
    i = int(row[0])
    j = int(row[1])
    if (i-8)%5==0 and (j-8)%5==0:
        points.append([i,j])
        
        kernel_Y_min = i - margin
        kernel_Y_max = i + margin + 1
        kernel_X_min = j - margin
        kernel_X_max = j + margin + 1
        
        discrepancy = int(row[-1])
        values.append(discrepancy)
        try:
            depth = (4.25) * (11 * 2.54 * 10) / (discrepancy*0.00349242)
        except:
            depth = np.inf

        for y in range(kernel_Y_min, kernel_Y_max):
            for x in range(kernel_X_min, kernel_X_max):
                discrepancies[y,x] = discrepancy
                depths[y,x] = depth

fig, ax = plt.subplots(ncols=2, figsize=(20,7))

ax[0].imshow(images[chosen_images[1]])
ax[0].set_title('Original Right Image')

cax = ax[1].imshow(discrepancies, cmap = 'rainbow')
ax[1].set_title('Disparity accros the image')
cbar = fig.colorbar(cax, ticks=np.arange(0,60,5))

fig, ax = plt.subplots(ncols=2, figsize=(20,7))

ax[0].imshow(images[chosen_images[1]])
ax[0].set_title('Original Right Image')

cax = ax[1].imshow(depths, cmap = 'rainbow')
ax[1].set_title('Depth accros the image')
cbar = fig.colorbar(cax, ticks=np.arange(0,100,10))

perlow, perhigh = np.percentile(np.nan_to_num(depths),[10,90])
depths = np.clip(depths, perlow, perhigh)

fig, ax = plt.subplots(ncols=2, figsize=(20,7))

ax[0].imshow(images[chosen_images[1]])
ax[0].set_title('Original Right Image')

cax = ax[1].imshow(depths, cmap = 'rainbow')
ax[1].set_title('Depth accros the image')
cbar = fig.colorbar(cax, ticks=np.arange(0,100))

from scipy.interpolate import griddata

points = np.array(points)
values = np.array(values)
grid_x, grid_y = np.mgrid[0:600, 0:800]
grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')

fig, ax = plt.subplots(ncols=2, figsize=(20,7))

ax[0].imshow(images[chosen_images[1]])
ax[0].set_title('Original Right Image')

cax = ax[1].imshow(grid_z0, cmap = 'rainbow')
ax[1].set_title('Interpolation of disparity accros the image')
cbar = fig.colorbar(cax, ticks=np.arange(-10,59,5))

val_per_low, val_per_high = np.percentile(values, [10,90])

depths_2 = []
for val in np.clip(values,val_per_low, val_per_high):
    depths_2.append((4.25 * 11 * 2.54 * 10) / (val * 0.00349242))
depths_2 = np.array(depths_2)
        
grid_x, grid_y = np.mgrid[0:600, 0:800]
grid_z0 = griddata(points, depths_2, (grid_x, grid_y), method='nearest')

fig, ax = plt.subplots(ncols=2, figsize=(20,7))

ax[0].imshow(images[chosen_images[1]])
ax[0].set_title('Original Right Image')

cax = ax[1].imshow(grid_z0, cmap = 'rainbow')
ax[1].set_title('Interpolation of disparity accros the image')
cbar = fig.colorbar(cax, ticks=np.arange(np.min(depths_2), np.max(depths_2), (np.max(depths_2)-np.min(depths_2))/10))

depths_3 = []
for val in values:
    depths_3.append((4.25 * 11 * 2.54 * 10) / (val * 0.00349242))

depth_arr = np.array(depths_3)
per_low, per_high = np.percentile(np.nan_to_num(depth_arr),[10,90])
depth_arr = np.clip(depth_arr, per_low, per_high)

grid_x, grid_y = np.mgrid[0:600, 0:800]
grid_z0 = griddata(points, depth_arr, (grid_x, grid_y), method='nearest')

fig, ax = plt.subplots(ncols=2, figsize=(20,7))

ax[0].imshow(images[chosen_images[1]])
ax[0].set_title('Original Right Image')

cax = ax[1].imshow(grid_z0, cmap = 'rainbow')
ax[1].set_title('Interpolation of disparity accros the image')
cbar = fig.colorbar(cax, ticks=np.arange(np.min(depth_arr), np.max(depth_arr), (np.max(depth_arr)-np.min(depth_arr))/10))

/Users/pthaper/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in double_scalars
  This is separate from the ipykernel package so we can avoid doing imports until